Simulating Highly Activated Sticking of H2 on Al(110): Quantum versus Quasi-Classical Dynamics

We evaluate the importance of quantum effects on the sticking of H2 on Al(110) for conditions that are close to those of molecular beam experiments that have been done on this system. Calculations with the quasi-classical trajectory (QCT) method and with quantum dynamics (QD) are performed using a model in which only motion in the six molecular degrees of freedom is allowed. The potential energy surface used has a minimum barrier height close to the value recently obtained with the quantum Monte Carlo method. Monte Carlo averaging over the initial rovibrational states allowed the QD calculations to be done with an order of magnitude smaller computational expense. The sticking probability curve computed with QD is shifted to lower energies relative to the QCT curve by 0.21 to 0.05 kcal/mol, with the highest shift obtained for the lowest incidence energy. Quantum effects are therefore expected to play a small role in calculations that would evaluate the accuracy of electronic structure methods for determining the minimum barrier height to dissociative chemisorption for H2 + Al(110) on the basis of the standard procedure for comparing results of theory with molecular beam experiments.

S2 S1. Setting up a 0 K metal slab for the convergence tests of the calculations with density functional theory.
All calculations with density functional theory (DFT) used for setting up the metal slab and for performing convergence tests for computing the molecule-metal surface interaction have been performed with the VASP 5.3.5 code 1-2 .
Calculations of H2 -Al(110) interaction energies are preceded by DFT calculations to set up the slab modeling the Al(110) surface according to the specific reaction parameter density functional (SRP DF) used. First, the bulk lattice constant al has been determined using a bulk lattice relaxation calculation. This calculation used a 1 x 1 x 1 super cell, a plane wave cutoff energy of 300 eV, 2 nd order Methfessel-Paxton smearing 3 with a smearing energy of 0.1 eV, and 25 x 25 x 25 Γ-centered k-points. Also, the calculations used projected augmented wave (PAW) pseudopotentials (PP) that were somewhat harder than the ordinary PAW PP available from VASP (i.e., we used the PP labeled Al_GW, with date stamp 19 March 2012, of VASP) as calculations with vdW-DF correlation functionals require somewhat harder PP than calculations with semi-local density functionals (DFs) within the generalized gradient approximation (GGA) would normally need [4][5][6][7] . Next, the distances between layers n and n+1 dn, n+1 (n=1 defines the surface layer) 10-layer Al slab (see below for the choice of 10 layers) in the slab were allowed to relax. These calculations employed the same PP as used to determine the bulk lattice constant, a 1 x 1 super cell, a plane wave cutoff energy of 300 eV, 2nd order Methfessel-Paxton smearing 3 with a smearing energy of 0.1 eV, 25 x 25 x 1 Γ-centered k-points, and a vacuum distance of 24 Å. The resulting "0 K" parameters are provided and compared with experiment in Table S1. Here, the experimental data for 0 K have been extracted from experiments presented in Ref. 8

as described
in the Supporting Information of Ref. 9 . The 0 K experimental data in Table S1 have not been   corrected for zero-point anharmonic expansion effects (they are the data from table S5 of Ref. 9 ,   S3 where the bulk lattice constant was obtained by multiplying the bulk interlayer distance determined in Ref. 8 S2. DFT convergence tests. The input parameters to the DFT calculations of the interaction energies were based on convergence tests for the BG1 and BG2 geometries of ( Fig.1) (i.e., the geometries labeled TS1 and TS2 in Ref. 9 ). These tests were performed using a slightly different functional as used in the calculation of the potential energy surface (PES), i.e., in Eq.1 68% RPBE 10 and 32% PBE 11 exchange was used in the convergence tests. The conclusions from the tests on BG2 were similar to those for BG1, therefore we only present results for BG1. The number of Al layers used in the Al slab to model the Al(110) surface was set to 10 in all tests, as previous convergence tests with the PBE DF had shown that DFT calculations with a 10-layer slab produced results that are very similar to the converged results obtained with a 19 layer slab 9 . For Al the calculations used the same PP as mentioned above for setting up the metal slab, and PAW PP for H that were somewhat harder than the ordinary PAW PP available from VASP (i.e., the PP labeled H_GW with date stamp 21 April 2008, of VASP), as calculations with vdW-DF correlation functionals require somewhat harder PP than calculations with semi-local GGA DFs would normally need [4][5][6][7] .
The calculations used a 3D Al lattice constant that was computed self-consistently with the DF used in the tests, and the same was true for the interlayer distances employed. The tolerance parameter governing the convergence of the DFT energy with electronic iteration number (EDIFF in VASP) was set to 10 -6 eV.
We first tested the convergence of the DFT results with the size of the surface unit cell (test series 1 in Table S2). Going from the (Nc x Nc) = (2x2) to the (3x3) surface unit cell changed the BG1 energy by no more than 0.2 kcal/mol. We chose the tighter (3x3) setting.

S4
In a second step (series 2 in table S2), we tested convergence with respect to the plane-wave cutoff EPW. Using EPW = 540 eV yields results that differ from results obtained with 700 eV by no more than 0.15 kcal/mol. We therefore adopted a value of 540 eV for EPW.
In the third step (series 3 in Table S2), we tested the number of k-points Nk required in the plane wave DFT calculations, using Monkhorst-Pack 12 k-point integration. (Is that correct, in the main paper it says we used a Γ-centered grid of k-points, is that the same?). Tests show that the BG1 energy changes by no more than 0.09 kcal/mol going from (Nk x Nk) = (8x8) k-points to (Nk x Nk) = (14x14) with the use of a (3x3) surface unit cell. We therefore chose to use a setup with (8x8) Monkhorst-Pack k-point integration.
We next (series 4 in Table S2) tested the convergence with respect to the vacuum distance DV.
Using DV = 16 Å yields results converged to within about 0.05 kcal/mol, and we therefore chose to use this value for DV.
Finally, we tested the convergence with respect to the energy smearing value Esm used with first order and second order Methfessel-Paxton smearing 3 (series 5 in Table S2). With either order the BG1 energy changes by less than 0.05 kcal/mol going from Esm = 0.2 to 0.1 eV. Also, the results for Esm = 0.1 eV are very similar for order 2 and order 1. We therefore chose to use Esm = 0.1 eV with first order Methfessel Paxton smearing.
The set up constituted by the input parameters to VASP is summarized in Table S3. S3. Setting up a 220 K slab for the computation of the potential energy surface.
To compute the DFT data on which the PES was based, we used a set up that is appropriate for the surface temperature Ts for which the dynamics calculations were performed, which was taken S5 equal to the Ts used in the experiments (220 K 13 ). Already in Ref. 9 , we determined the experimental bulk interlayer distance at 0 K by extrapolating the measured temperature dependent bulk interlayer distance to 0 K, yielding a value of 1.4239 Å (we used the data in fig.6 of Ref. 8 for this, see also Table S1). Multiplying this value with 2 2 yields a 0 K bulk lattice constant of 4.0274 Å. Next, we determined the slope of the temperature dependent bulk interlayer distance from the experimental data, yielding 4.3 x 10 -5 Å K -1 (again, we used the data in fig.6 of Ref. 8 for this). Using these two data yields a value of 1.4335 Å for the experimental bulk interlayer distance at 220 K. Multiplying this value with 2 2 yields an experimental 220 K bulk lattice constant of 4.0544 Å.
We then obtain the 220 K "DFT lattice constant" using Substituting the 0 K lattice constant obtained from DFT (of 4.08708 Å, see Table S1), we obtain a value of 4.11448 Å for a l DFT (220 K) . In an entirely analogous fashion, a 220 K DFT bulk inter layer distance d n,n+1 DFT (220 K) of 1.4547 Å is obtained.
We next obtained interlayer distances between layers 1 and 2, and between layers 2 and 3, from our DFT calculations and the experiments of Ref. 8 . First we determined experimental 0 K interlayer distances using extrapolation of the temperature dependent interlayer distances displayed in fig.6 of Ref. 8  and d 5,6 change with temperature. We therefore assume that these can be obtained from the bulk interlayer distance DFT value d n,n+1 DFT (220 K) , the percentage changes of these interlayer distances in the relaxed 0 K DFT slab, and d n,n+1 DFT (0 K) . The percentage changes of these 0 K interlayer distances, which are tabulated in Table S1, were Δ(d 3,4 ) = -3.2%, Δ(d 4,5 ) = 1.9%, and Δ(d 5,6 ) = 0.8, respectively. The resulting 220 K values of the corresponding interlayer distances are presented in Table S1.
S4. Input parameters to calculations with the time-dependent wave packet method.
Calculations with the time-dependent wave packet (TDWP) method were performed for incidence energies between 0.05 and 1.05 eV. It is not possible to obtain results for this entire energy range in one calculation: it would not be possible to contain the initial wave packet in the appropriate range of momenta directed at the surface and to obtain accurate results at the same time. Therefore, this energy range was split into four smaller energy ranges, and TDWP calculations were performed for these ranges (see Table S4) separately. In each calculation, the Z-dependence (r-dependence) of the wave function was described with a Fourier grid starting at the start value Z start ( r start ) provided in Table S4 with N Z ( N r ) subsequent grid points with grid spacings of ΔZ ( Δr ), as all indicated in Table S4 (the coordinates used have been specified in Fig.1 of the main paper). A projection operator formalism 14 was used to bring in the initial wave packet on a separate, long one-dimensional grid in order to be able to reduce the grid size in Z that is associated with the large scattering basis set. In each calculation this grid has the same start value and grid spacings as the ordinary Z-grid, but it has a larger number of points N Z sp (see Table S4). Furthermore the scattering basis set has periodic Fourier grids in X and Y with N X and N Y points, respectively, and in the finite-basis representation (FBR) spherical harmonics are S7 employed with the maximum value of j equal to j max and the maximum value of m j equal to m j max as also given in Table S4.
On the grids in Z and r defined above complex absorbing potentials (CAPS) are defined in Table   4 from a starting value to an end value, and these are stated for the regular Z-grid, the Z-grid used to bring in the initial wave function with the projection operator formalism, and the r-grid using parameters with names that are self-evident. The strength of the squared negative imaginary potentials used as CAPs was taken in such a way that optimal absorption takes place for the translational energy indicated in Table S4 (in eV), using the theory of Ref. 15 and also based on the range of the grid over which the optical potential is defined.
The wave function was propagated with a time step Δt up to a final time t f selected in such a way that the remaining norm of the wave function on the grid was in all cases less than 0.2%.
The initial wave function was centered on the value of Z 0 indicated. Once again these parameters are provided in Table S4. S8 Table S1. Parameters characterizing the Al(110) slab obtained with the SRP DF in calculations using the static surface approximation and in calculations used to set up a potential energy surface for H2 + Al(110) for a surface temperature of 220 K. Experimental parameters (in brackets) are taken from the experimental data presented in Ref. 8 , and from the analysis performed of these data in Ref. 9 . All interlayer distances are in Å.  Figure S1. Sticking probabilities computed with the QCT method using the NMC (blue squares) and FMC (green triangles) procedures, as described in the text of the main paper, are compared. Incidence energy [kcal/mol]